Setup
# Remove all from the current workspace
rm(list = ls())
# Set the working directory
setwd('C:\\Users\\srica\\Desktop\\Radio\\projekt\\PDR Tema 9')
# Create a directory named 'plots' if it doesn't already exist
if (!file.exists('plots')) dir.create('plots')
# Create a directory named 'results' if it doesn't already exist
if (!file.exists('results')) dir.create('results')
# Import the 'stringr' library
library(stringr)
Constants
# Set the necessary constants
constants <- data.frame(
fiPole = 78.3 * pi / 180,
lambdaPole = 291 * pi / 180,
ippHeight = 350,
earthRadius = 6378,
lightSpeed = 2.99792458e+08
)
print(constants)
## fiPole lambdaPole ippHeight earthRadius lightSpeed
## 1 1.366593 5.078908 350 6378 299792458
User coordinates
# receiver (user) cooridinates (Pevex Kukuljanovo)
userFiLambda <- data.frame(
fiUser = 45.331360488178596*(pi/180), # longitude
lambdaUser = 14.509391791906435*(pi/180) # latitude
)
print(userFiLambda)
## fiUser lambdaUser
## 1 0.7911815 0.2532367
# real Cartesian receiver (user) coordinates
userCoords <- data.frame(
X = constants$earthRadius * cos(userFiLambda$lambdaUser) * cos(userFiLambda$fiUser),
Y = constants$earthRadius * cos(userFiLambda$lambdaUser) * sin(userFiLambda$fiUser),
Z = constants$earthRadius * sin(userFiLambda$lambdaUser)
)
print(userCoords)
## X Y Z
## 1 4340.767 4391.267 1597.936
Alpha and beta
# Klobuchar ionospheric model constants
ionisation <- data.frame(
alphaIon <- c(1.4900e-08, -7.4510e-09, -5.9600e-08, 1.1920e-07),
betaIon <- c(1.2900e+05, -1.9660e+05, 6.5540e+04, 3.2770e+05)
)
print(ionisation$alphaIon)
## [1] 1.490e-08 -7.451e-09 -5.960e-08 1.192e-07
print(ionisation$betaIon)
## [1] 129000 -196600 65540 327700
sp3 data
# Defining variables used in Klobuchar model and loading data
firstDate <- NA
sp3DataFormatted <- list()
sp3Data <- read.delim('data\\data.SP3', header=FALSE)
# Looping over data (skipping first 23 lines which are header data) and loading it into an array
for (rowIndex in 23:nrow(sp3Data)-1) {
# Loading data from sp3Data at index 'rowIndex' and splitting the contents (after removing whitespaces) into an array
# gsub replaces all types of whitespaces with ' '
# str_trim removes leading and trailing whitespaces
# str_split splits the content into an array depending on ' '
row <- str_split(str_trim(gsub('\\s+', ' ', sp3Data[rowIndex, ])), ' ')[[1]]
# If row length is 7 that line contains timestamp for the next 78 observations
# Timestamp format: * year month day hour minutes seconds
if (length(row) == 7) {
# Load timestamp into a variable
dateTime <- strptime(paste(row[2], row[3], row[4], row[5], row[6]), format = "%Y %m %d %H %M")
# Save first observation into a separate variable for later calculations
if (is.na(firstDate)) firstDate <- dateTime
# gpsTime variable contains difference between first reading and current reading
gpsTime <- as.numeric(difftime(dateTime, firstDate, units='mins'))
# Skips to the next step of the for loop
next
}
# All other rows contain observations for each satellite
# Observation format: id x y z clock deviations
data = data.frame(DateTime=dateTime,
GPSTime=gpsTime,
X=as.numeric(row[2]),
Y=as.numeric(row[3]),
Z=as.numeric(row[4]))
# Save observation into an array
sp3DataFormatted[[row[1]]] <- rbind(sp3DataFormatted[[row[1]]], data)
}
# Cleaning unneeded variables from the current workspace
rm(rowIndex, row, data, dateTime, gpsTime, firstDate)
Data calculations and saving
# Calculate Klobuchar ionospheric delay
getIonDist <- function(gpsData){
# Difference between the user and satellite height
elevation <- gpsData$Z - userCoords$Z
# Difference between the satellite and user coordinates
# atan() function is used to compute arctangent angle
azimuth <- atan((gpsData$X - userCoords$X) - (gpsData$Y - userCoords$Y))
# Angle between the Earth and satellite
psi <- pi / 2 - elevation - asin(constants$earthRadius / (constants$earthRadius + constants$ippHeight) * cos(elevation))
# Determine the latitude of the IPP point
fiIon <- asin(sin(userFiLambda$fiUser) * cos(psi) + cos(userFiLambda$fiUser) * sin(psi) * cos(azimuth))
# Determine the longitude of the IPP point
lambdaIon <- userFiLambda$lambdaUser + psi * sin(azimuth) / cos(fiIon)
# Determine the geomagnetic latitude of the IPP point
fiMag <- asin(sin(fiIon) * sin(constants$fiPole) + cos(fiIon) * cos(constants$fiPole) * cos(lambdaIon - constants$lambdaPole))
# Determine the local time in the IPP point
time <- 43200 * lambdaIon / pi + gpsData$GPSTime
if (time >= 86400) time <- time - 86400
if (time < 0) time <- time + 86400
# Calculate the amplitude of the ionospheric delay
azimuthIon <- 0
for (i in seq(1, 4, 1)) azimuthIon <- azimuthIon + ionisation$alphaIon[i] * `^`(fiMag / pi, i - 1)
if (azimuthIon < 0) azimuthIon <- 0
# Calculate the period of the ionospheric delay
psiIon <- 0
for (i in seq(1, 4, 1)) psiIon <- psiIon + as.double(ionisation$betaIon[i]) * `^`(fiMag / pi, i - 1)
if (psiIon > 72000) psiIon <- 72000
# Calculate the phase of the ionospheric delay
xIon <- 2 * pi * (time - 50400) / psiIon
# Calculate ionospheric mapping function (inclination factor)
Fun <- `^`(1 - `^`(constants$earthRadius / (constants$earthRadius + constants$ippHeight) * cos(elevation), 2), -1 / 2)
# Calculate the value of vertical ionospheric delay
if (abs(xIon) < pi / 2) dIon <- (5e-9 + azimuthIon * cos(xIon)) * Fun
else dIon <- 5e-9 * Fun
# Calculate the distance from ionospheric delay by multiplying ionospheric delay with the speed of light
dIon <- dIon * constants$lightSpeed
# Calculate the Euclidean distance in km from the actual distance between the user and the GPS satellite
diffXSquared <- (userCoords$X - gpsData$X) ** 2
diffYSquared <- (userCoords$Y - gpsData$Y) ** 2
diffZSquared <- (userCoords$Z - gpsData$Z) ** 2
euclidianDistance <- sqrt(diffXSquared + diffYSquared + diffZSquared)
# Calculate the pseudo distance(m) from the Euclidean distance(km) and ionospheric distance(m)
distance <- euclidianDistance * 1000 + dIon
return(list(dIon=dIon, distance=distance))
}
# Function for saving and plotting data
saveData <- function(gpsName, lines, ionVector, distVector){
time <- seq(0, 24, len=length(ionVector))
# Create files for the current GPS
fileIon <- file(paste('results/', gpsName, 'ionOutput.txt'))
fileDist <- file(paste('results/', gpsName, 'distOutput.txt'))
# Save data into the file
writeLines(as.character(ionVector), fileIon)
close(fileIon)
writeLines(lines, fileDist)
close(fileDist)
# Plotting data
plot(time, ionVector, type='l', main=str_interp(paste('Klobucharov model za', gpsName)), xlab='time[h]', ylab='dIon[m]')
plot(time, distVector, type='l', main=str_interp(paste('Pseudoudaljenost za', gpsName)), xlab='time[h]', ylab='distance[m]')
# Saving plots
png(paste('plots/', gpsName, 'ionPlot.png'), width=20, height=10, units='cm', res=500)
plot(time, ionVector, type='l', main=str_interp(paste('Klobucharov model za', gpsName)), xlab='time[h]', ylab='dIon[m]')
dev.off()
png(paste('plots/', gpsName, 'distPlot.png'), width=20, height=10, units='cm', res=500)
plot(time, distVector, type='l', main=str_interp(paste('Pseudoudaljenost za', gpsName)), xlab='time[h]', ylab='distance[m]')
dev.off()
}
# Preparing data to save
emptySpaces <- strrep(' ', 25)
timestampHeader <- substr(str_interp('TIME [min]${emptySpaces}'), 0, 15)
distanceHeader <- substr(str_interp('DISTANCE [m]${emptySpaces}'), 0, 15)
idHeader <- substr(str_interp('ID ${emptySpaces}'), 0, 10)
header <- str_interp('${timestampHeader} ${idHeader} ${distanceHeader}')
# Looping over each satellite
for(gpsName in names(sp3DataFormatted)){
lines <- header
ionVector <- c()
distVector <- c()
# Looping over data for each satellite and formatting the output
for(gpsIndex in 1:nrow(sp3DataFormatted[[gpsName]])){
gpsData <- sp3DataFormatted[[gpsName]][gpsIndex, ]
ionDist <- getIonDist(gpsData)
ionVector <- append(ionVector, ionDist$dIon)
distVector <- append(distVector, ionDist$distance)
gpsTimeFormatted <- substr(str_interp('${gpsData$GPSTime}${emptySpaces}'), 0, 15)
distanceFormatted <- substr(str_interp('${ionDist$distance}${emptySpaces}'), 0, 15)
gpsNameFormatted <- substr(str_interp('${gpsName}${emptySpaces}'), 0, 10)
lines <- append(lines, str_interp('${gpsTimeFormatted} ${gpsNameFormatted} ${distanceFormatted}'))
}
# Saving data using the previously defined function
saveData(gpsName, lines, ionVector, distVector)
print(paste('Finished', gpsName))
}


## [1] "Finished PG01"


## [1] "Finished PG02"


## [1] "Finished PG03"


## [1] "Finished PG04"


## [1] "Finished PG05"


## [1] "Finished PG06"


## [1] "Finished PG07"


## [1] "Finished PG08"


## [1] "Finished PG09"


## [1] "Finished PG10"


## [1] "Finished PG11"


## [1] "Finished PG12"


## [1] "Finished PG13"


## [1] "Finished PG14"


## [1] "Finished PG15"


## [1] "Finished PG16"


## [1] "Finished PG17"


## [1] "Finished PG18"


## [1] "Finished PG19"


## [1] "Finished PG20"


## [1] "Finished PG21"


## [1] "Finished PG22"


## [1] "Finished PG23"


## [1] "Finished PG24"


## [1] "Finished PG25"


## [1] "Finished PG26"


## [1] "Finished PG27"


## [1] "Finished PG28"


## [1] "Finished PG29"


## [1] "Finished PG30"


## [1] "Finished PG31"


## [1] "Finished PG32"


## [1] "Finished PR01"


## [1] "Finished PR02"


## [1] "Finished PR03"


## [1] "Finished PR04"


## [1] "Finished PR05"


## [1] "Finished PR07"


## [1] "Finished PR08"


## [1] "Finished PR09"


## [1] "Finished PR11"


## [1] "Finished PR12"


## [1] "Finished PR13"


## [1] "Finished PR14"


## [1] "Finished PR15"


## [1] "Finished PR16"


## [1] "Finished PR17"


## [1] "Finished PR18"


## [1] "Finished PR19"


## [1] "Finished PR20"


## [1] "Finished PR21"


## [1] "Finished PR24"


## [1] "Finished PR25"


## [1] "Finished PE02"


## [1] "Finished PE03"


## [1] "Finished PE04"


## [1] "Finished PE05"


## [1] "Finished PE07"


## [1] "Finished PE08"


## [1] "Finished PE09"


## [1] "Finished PE10"


## [1] "Finished PE11"


## [1] "Finished PE12"


## [1] "Finished PE13"


## [1] "Finished PE14"


## [1] "Finished PE15"


## [1] "Finished PE18"


## [1] "Finished PE19"


## [1] "Finished PE21"


## [1] "Finished PE24"


## [1] "Finished PE25"


## [1] "Finished PE26"


## [1] "Finished PE27"


## [1] "Finished PE30"


## [1] "Finished PE31"


## [1] "Finished PE33"


## [1] "Finished PE34"


## [1] "Finished PE36"
# Cleaning variables from the current workspace
rm(distanceFormatted, distanceHeader, emptySpaces, gpsIndex, gpsName, gpsTimeFormatted, gpsNameFormatted, idHeader, header, ionVector, distVector, lines, timestampHeader)